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ABSTRACT 

We perform three-dimensional cosmological simulations to examine the growth of 
metal-free, Population III (Pop III) stars under radiative feedback. We begin our 
simulation at z = 100 and trace the evolution of gas and dark matter until the forma- 
tion of the first minihalo. We then follow the collapse of the gas within the minihalo 
up to densities of n — 10 12 cm' 3 , at which point we replace the high-density parti- 
cles with a sink particle to represent the growing protostar. We model the effect of 
Lyman- Werner (LW) radiation emitted by the protostar, and employ a ray-tracing 
scheme to follow the growth of the surrounding H n region over the next 5000 yr. We 
find that a disk assembles around the first protostar, and that radiative feedback will 
not prevent further fragmentation of the disk to form multiple Pop III stars. Ionization 
of neutral hydrogen and photodissociation of H2 by LW radiation leads to heating of 
the dense gas to several thousand Kelvin, and this warm region expands outward at 
the gas sound speed. Once the extent of this warm region becomes equivalent to the 
size of the disk, the disk mass declines while the accretion rate onto the protostars 
is reduced by an order of magnitude. This occurs when the largest sink has grown 
to ~ 20 M while the second sink has grown to ~ 7 M , and we estimate the main 
sink will approach an asymptotic value of 30 M Q by the time it reaches the main 
sequence. Our simulation thus indicates that the most likely outcome is a massive 
Pop III binary. However, we simulate only one minihalo, and the statistical variation 
between minihaloes may be substantial. If Pop III stars were typically unable to grow 
to more than a few tens of solar masses, this would have important consequences for 
the occurence of pair-instability supernovae in the early Universe as well as the Pop 
III chemical signature in the oldest stars observable today. 

Key words: stars: formation - Population III - galaxies: formation - cosmology: 
theory - first stars - early Universe 



1 INTRODUCTION 

The first stars were the earliest luminous objects to form at 
the end of the 'Dark Ages' that followed the emission of the 
Cosmic Microwave Background. These stars are thought to 
have formed around z > 20 within minihaloes of mass M ~ 
10 6 Mq , when the dark matter (DM) potential well of the 
mini halo becomes large enough to gather in primordial gas 
(e.g.lHaiman et al.ll99r3 : ITegmark et al.lll997l ; lYoshida et all 



Kitayama et all |2004 ISokasian et all |2004 IWhalen et all 
2004 1 Alvarez et al.ll2006l ; I Johnson et al.ll2007l ). With the re- 



lease of the first metals through their possible supernova 
(SN) deaths, they also provided the early metal e nrich- 



12003 ). Also known as Population III (Pop I II) stars, they are 
the early drivers of cosmic evolution (e.g. [B arkana fc Loeb 
200l| : lBromm fc Larsonll2004ICiardi fc Ferrarall2005l ; lGlovei 
20051 : iBromm et al.ll2009l : lLoebll2010l ). Through their emis- 
sion of ionizing radiation over their lifetime, Pop III stars 
are responsible for the beginning of cosmic reionization (e.g. 
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Mori et all [200I 1 Bromm et all 120031 IWada & Venkatesan 


2003: Norman et al.ll2004 


;lTornatore et al|2007t ICreif et al. 


2001 boirt Wise & Abell 


200S; Maio et alKOllI; recently re- 


viewed in Karlsson et al. 


20111). 



The extent to which Pop III stars can modify their sur- 
roundings is crucially dependent upon their mass, as this 
is the main characteristic that determines the star's lumi- 
nosity and ionizing radiation output. Furthermore, their 
mass determines the typ e of stellar death they will un- 



■yp e 

dergo (|Heger et al.l l2003h . Only stars in the mass range 
of 140 Mq < M* < 260 M are predicted to explode a s 
pair-instability supernovae (PISNe uHeger fc Wooslevll2002l ). 
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while stars with masses in the range 40 Mq < M* < 140 Mq 
are thought to collapse directly into black holes. Below 40 
Mq, stars are again expected to explode as core-collapse 
SNe, leaving behind a neutron star or black hole. 

Previous work has found that Pop III stars begin 
as v ery small protostars of initial mass ~ 5 x 10 -3 
M (|Omukai fc Nishil Il998l ; lYoshida et all |200S|) . Contin- 
ued accretion onto the protostar over time ultimately leads 
to stars significantly larger than the ini tial seeds (e.g. 
lOmukai fc Pallall2003l ;lB romm & Locb 2004). The final mass 
reached by Pop III stars therefore depends on the rate and 
duration of the accretion. Early numerical studies found that 
Pop III stars are likely t o reach very high m asses (> 100 
M ; e.g. lAbel et~ai1l2002l : iBromm et alJl2002l ). The lack of 
metal and dust cooling in primordial gas leads to higher tem- 
peratures and greater accretion rates as compared to cur- 
rent star-forming regions such as the giant molecular clouds 
within the Milky Way. 

The accretion history of Pop III stars depends critically 
on the radiative feedback exerted during their growth phase, 
an effect not included in the earliest three-dimensional simu- 
lations. The strength of protostellar feedback in turn hinges 
upon the thr ee-dimensional structure of the surrounding gas. 
For instance, lOmukai fc Inutsukal (|2002T l find that for spher- 
ically symmetric accretion onto an ionizing star, the forma- 
tion of an H II region in fact does not impose an upper 
m ass limit to Pop III st ars. Similarly, the analytical study 
bv lMcKee fcTa n (2008) found that Pop III stars can grow 
to greater than 100 Mq even as a protostar's own radiation 
ionizes its surroundings. Although the ionization front will 
expand along the polar regions perpendicular to the disk, 
mass from the disk itself can continue to accrete onto the 
star until this is halted through photoevaporation. Recent 
studies of present-day star formation also support the pic- 
ture that stars can reach very high masses through disk ac- 
cretio n. Fo r instance, numerica l studies by iKrumholz et all 
(2009) and lKuiper et al.l (|201ll ) both find that strong feed- 
back from radiation pressure will not halt mass flow through 
the stellar disk before the star has reached > 10 Mq. The 
final stellar masses are likely even greater, though the sim- 
ulations were not followed for sufficiently long to determine 
this. There is also growing observational evidence that mas- 
sive star-forming regions exhibit disk structure and rota- 
tional motion ( e.g. ICesaroni et al ] |2007l ;l Beuther et al . 2009; 
iLiu et alJlioTl . 

Meanwhile, the picture of a single massive Pop III star 
forming in a minihalo has been complicated by more re- 
cent work. Simulations by Clark et al. (2008, 2011a) em- 
ploying idealized initial conditions found that primordial 
star-forming gas can undergo fragmentation t o form Pop 
III m ultipl e systems, while t he simulations of I Turk et al.l 
(2009) and IStacv et al.l (|2010h established such fragmenta- 
tion also when initialized on cosmological scales. Further 
work revealed that gas fragmentation can occur even on 
very small scales (~ 10 AU) and in the majority of min i- 
haloes, if not nearly all (Clark et al. 2011b. lGreif et al.ll201l[) . 
These studies tentatively imply that the typ ical Pop HI 
mass may be somewhat lower than ~ 100 Mq . ISmith et al.l 
Hon]) recently found that, even under feedback from pro- 
tostellar accretion luminosity, such fragmentatio n may be 
reduced but not halted. While [Smith et all (|201 lh included 
a heating term derived from the accretion luminosity, they 



did not explicitly account for molecular photodissocation 
by Lyman- Werner (LW) radiation from the protostar, and 
they did not include the effects of ionizing radiation that 
will become important once the stars have grown to larger 
m asses (> 10 Mq). The r ecent two-dimensional calculation 
bv lHosokawa et all (|201 lh found that Pop III stars will grow 
to 40 Mq, after which accretion will be shut off by radia- 
tive feedback. However, the effects of three-dimensional non- 
axisymmetry could not be addressed in this study. 

Whether Pop III stars can attain very large masses un- 
der feedback and while within a multiple system is pivotal 
to understanding their potential for large-scale feedback ef- 
fects, such as the suppresion or enhancement of the star- 
formation rate in neighboring minihaloes. It also determines 
whether they may be observed as ga mma-ray bursts (GRBs) 
or extremely energetic PISNe (e.g. IBromm fc Loebl |2002|. 
| 2006t iGou et all [20041 ; iBelczvnski et all 120071 ; IStacv et al l 
120111 ). Furthermore, the fragmentation seen in recent work, 
along with the possible ejection of low- mass Pop HI stars 
from their host s tar-forming disks (e.g. iGreif et al.l l201ll ; 
ISmith et ai]|201ll ), opens the possibility that small, long- 
lived Pop III stars may still be observed today. This de- 
pends, however, on uncertain factors such as the final masses 
reached by the ejected stars and the amount of metal- 
enriched material accreted at later times while being incor- 
porated int o larger galaxies, which could mask the stars as 
Pop II fe.g. lFrebel et al.ll2009l ; [Johnson fc Khochfaijl201lT ). 

To further explore the range of masses possible for Pop 
III stars, we perform a three-dimensional cosmological sim- 
ulation to study the feedback effects of a protostar on its 
own accretion and on further fragmentation within its host 
minihalo. We initialize the simulation with sufficient reso- 
lution to follow the evolution of the star-forming gas up to 
densities of 10 12 cm -3 . At this density we employ the sink 
particle method, allowing us to study the subsequent disk 
formation and fragmentation of the gas over the following ~ 
5000 yr. We include H2 dissociating LW feedback from the 
most massive star in the simulation, and we use a ray-tracing 
scheme to follow the growth of the star's H 11 region once it 
has become massive enough to ionize the surrounding gas. 
We compare this to a simulation with the same initialization 
but no radiative feedback. This allows for a direct evaluation 
of how radiative feedback alters the mass growth of Pop III 
stars, the rate of disk fragmentation, and the formation of 
additional stars within the disk. We describe our numerical 
methodology in Section 2, while in Section 3 we present our 
results. We conclude in Section 4. 



2 NUMERICAL METHODOLOGY 
2.1 Initial Setup 

We ran our simulation using GADGET 2, a widely-teste d 
three-dimensional N-body and SPH code (|Springej|2005h . 
We initialized our si mulation using a s napshot from the 
previous simulation of IStacv et al.l (|2010T ). In particular, we 
chose the snapshot in which the dense gas in the center of 
the minihalo has first reached 10 s cm -3 . This simulation 
was originally initialized at z = 100 in a periodic box of 
length 140 kpc (comoving) using both SPH and DM parti- 
cles. This was done in accordance with a ACDM cosmology 
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with ft A = 0.7, fi M = 0.3, n B = 0.04, and H = 70 km 
s _1 Mpc -1 . To accelerate structure formation, we used an 
artificially enhanced no rmalization of the power spectrum, 
erg = 1.4. As discussed in lStacv et al.l |20ld ). we verified that 
the density and velocity fields in the center of the minihalo 
closely resembled those in previous simulations. We further- 
more found that the angular momentum profile of the mini- 
halo gas immediately before initial sink f ormation w a s very 
simil ar to the cosmologies ! simulations of lAbel et all l|2002n 
and lYoshida et all (^OOd - ). despite their lower values of ae 
(0.7 and 0.9, respectively). 

The high resolution of our simulation was achieved 
through a stand ard hierarchical zoom- in procedure (see 
IStacv et all l20ld for more details). We added three addi- 
tional nested refinement levels of 40, 30, and 20 kpc (comov- 
ing), centered on the location where the first minihalo will 
form. At each level of refinement, we replaced every 'parent' 
particle with eight 'child' particles, such that at the high- 
est refinement level each parent particle has been replaced 
with 512 child particles. The highest-resolution gas particles 
have a mass jtisph = 0.015 M©, so that the mass resolution 
is M les ~ l.SA^neighmspH < 1 Mq, where A ne igh — 32 is the 
typic al number of particle s in the SPH smoothing kernel 
(e.g. iBate fc Burkertlll997f l. 

2.2 Cut-Out Technique 

Inclusion of feedback significantly reduces the simulation 
timesteps as compared to our control case with no feed- 
back. To facilitate faster computation speeds in our 'with- 
feedback' case, once the main sink becomes massive enough 
to emit ionizing radiation we implement a 'cut-out' tech- 
nique. In particular, we remove from the simulation box all 
particles that are located beyond 10 pc (physical) from the 
main sink, thereby following only the gravitationally bound 
central gas and discarding the very slowly evolving outer re- 
gions. The cut-out region corresponds to the central < 4000 
Mq of gas. This technique reduces total computation time 
by nearly an order of magnitude while having a minimal 
effect on the simulated accretion history of the main sink, 
since by this point the central gas is dense, self-gravitating 
and no longer influenced by the gravity of the outer minihalo 
or DM on larger scales. 

We furthermore note that the gas at the edge of the 10 
pc 'cut-out' region has a typical density of ~ 10 2 cm -3 , 
corresponding to a free-fall time of ~ 10 7 yr, and thus 
should undergo little evolution within our simulation time of 
5000 yr. In addition, though our method leads to a vacuum 
boundary condition at the edge of the cut-out, this should 
not be problematic. The boundary conditions may lead to 
the propagation of a rarefaction wave starting from the cut- 
out edge, but this will only travel a distance of c s t, where 
c s is the gas soundspeed (~ 2 km s _1 ), and the time t is 
5000 yr. This corresponds to a distance of ~ 10 -2 pc (2000 
AU) from the cut-out edge, a very small distance compared 
to the 10 pc box size. 

2.3 Chemistry, Heating, and Cooling 

W e use the same che mistry and cooling network as described 
ni iGreif et all j2009l V The code follows the abundance evo- 
lution of H, H+, FT, H 2 , Uj, He, He+, He++, and e~, as 



well as the three deuterium species D, D + , and HD. All 
relevant cooling mechanisms are accounted for, including 
H2 cooling through collisions with He and H atoms and 
other H2 molecules. Also included are cooling through H 
and He collisional excitation and ionization, recombination, 
bremsstrahlung, and inverse Compton scattering. We finally 
note that H2 cooling through collisions with protons and 
electrons is i ncluded, as they play an important role within 
H 11 regions (|Glover fc Abell2008h . 

One important difference we note between the high den- 
shy in > 10 9 cm" 3 ) evolutio n in the current simulation 
and that of IStacv et all l|2010h is that we here account for 
the optical thickness of the H2 ro-vibrational lines, which 
reduces the effectiveness of these lines in cooling the gas. 
We include this effect using the Sobol ev approximation (see 
lYoshida et aUl200d ; IGreif et al.ll201ll for more details). For 
the three-body reactions 

H + H + H^Ha + H 

and 

H + H + H 2 ->■ H 2 + H 2 

we ch oose the rate coefficients adopted by iPalla et al.l 
|l983t ). We note, however, that these re action rates are s till 
subject to significant uncertainties fsee lTurk et aljfeoill for 
a discussion). 

2.4 Sink Particle Method 

When an SPH particle reaches a density of n max = 10 12 
cm" 3 , we convert it to a sink particle. If a gas particle is 
within the accretion radius r acc of the sink, and if it is not 
rotationally supported against infall onto the sink, the sink 
accretes the particle. The sink thus accretes the particles 
within its smoothing kernel immediately after it first forms. 
We check for rotational support by comparing the angular 
momentum of the nearby gas particle, j'sph = v TO td, with 
the angular momentum required for centrifugal support, 
Jccnt = V GA/ S inkr a cc , where v ro t and d are the rotational 
velocity and distance of the particle relative to the sink. If 
a gas particle satisfies both d < r acc and jspu < jcont , it is 
removed from the simulation, and the mass of the accreted 
particle is added to that of the sink. 

Our sink accretion algorithm furthermore allows for the 
merging of two sink particles. We use similar criteria for sink 
particle merging as for sink accretion. If the smaller, sec- 
ondary sink is within r acc of the more massive sink and has 
specific angular momentum less than j ce nt of the larger sink, 
the sinks are merged. After an accretion or merger event, the 
position of the sink is set to the mass-weighted average of 
the sink's former position and that of the accreted gas or 
secondary sink. The same is done for the si nk velocity. We 
note that, as discussed in IGreif et alj (|201ll ). modifications 
to the sink merging algorithm can significantly alter the sink 
accretion history. Future work will include studies of how a 
different technique for sink merging would modify our re- 
sults. 

We set the accretion radius equal to the resolution 
length of the simulation, r acc = L rcs ~ 50 AU, where 

L TeB ~ 0.5 1 , (1) 
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with p m ax — limaxfflH and tuh being the proton mass. The 
sink particle's mass, M s i n k, is initially close to the resolu- 
tion mass of the simulation, M Tes ~ 0.7 Mg. Sink particles 
are held at a constant density and temperature of n max = 
10 12 cm' 3 and 650 K, such that their pressure is also kept to 
the corresponding value. Providing the sink with a tempera- 
ture and pressure prevents the existence of a pressure deficit 
around the sink, which wo uld otherwise lea d to artificially 
high accretion rates (e .g. iBate et"al] 1 19951 : iBromm et al.l 
120021 ; iMartel et al.ll2006T ). The sink particles do continue to 
evolve in position and velocity, however, through gravita- 
tional and hydrodyn amic interactions. 

As discussed in IBromm et alj (|2002T l and IStacv et alj 

(2010), our sink formation criteria well represent regions 
that will truly collapse to stellar densities. Before cross- 
ing the density threshold to become a sink, a gas particle 
must collapse two orders of magnitude above the average 
disk density, ~ 10 10 cm" 3 . Our high density threshold and 
small value for r acc ensure that sinks are formed only from 
gravitationally collapsing gas. 

Following the long-term evolution of the star-forming 
gas would not be feasible without the sink particle method. 
By preventing gas evolution to ever higher densities, we 
avoid the problem of increasingly small numerical timesteps, 
a problem also known as 'Courant myopia.' We thus are 
able to see how the surrounding region of interest evolves 
over many dynamical times. With sink particles, we can 
furthermore bypass the need to incorporate the chemistry, 
hydrodynamics and radiative transfer that comes into play 
at extremely high densities (n > 10 12 cm -3 ). Finally, sink 
particles provide a convenient way to directly measure the 
accretion rate onto the protostellar region. 



2.5 Ray-tracing Scheme 

Once the first sink is formed, this particle is used as the 
source of protostellar LW and ionizing radiation. While 
the protostar is less massive than 10 M©, LW radiation is 
the only source of feedback. After the protostar is massive 
enough to emit ionizing radiation, however, a compact H 11 
region develops. We model the growth of the surrounding 
I- front using a ray-tr acing scheme which closely follows that 
of iGreif et ail £2009). A spherical grid with ~ 10 5 rays and 
200 radial bins is then created around the sink particle. The 
minimum radius is determined by the distance between the 
sink and its closest neighboring SPH particle, and we up- 
date this structure each time the ray-tracing is performed. 
Because the sink accretes most particles within r acc = 50 
AU, the minimum radius is usually > 50 AU. The maxi- 
mum radius is chosen as 10 pc (physical), the size of the 
cut-out simulation. This value easily encompasses the entire 
H 11 region in our simulation. The radial bins within 75 AU 
of the minimum radius are spaced at intervals of 1.5 AU. 
Outside this distance the bins are logarithmically spaced. 
The location of each particle is then mapped onto the cor- 
responding bin within the spherical grid, and each particle 
contributes its density and chemical abundances to the bin 
proportional to its density squared. 

Next, the ionization front equation is solved along each 

ray: 



n n r t 



dt 



4-7T 



ob / n e n + r dr 



(2) 



where n n , n e , and n+ are the number densities of neutral 
particles, electrons, and positively charged ions, respectively. 
The location of the ionization front with respect to the sink 
is denoted by n . iVion is the number of ionizing photons emit- 
ted per second, and as is the case B recombination coeffi- 
cient. We use as = 1.3 x 10 -12 cm 3 s _1 for He ill recombina- 
tions to He 11, and «b = 2.6x 10" 13 c m 3 s" 1 for He 11 and H 11 
recom binations to the ground state (|Osterbrock fc Ferlandl 
l2006h . 

The above I-front equation assumes that the ioniza- 
tion front is expanding into neutral, non-molecular gas. This 
turns out to be an accurate assumption despite that the 
star is surrounded by a molecular disk. As will be further 
described below, the ionization front expands only into the 
lower-density polar regions which are indeed non-molecular. 
Note that the gas does not become fully molecular until 
it reaches densities of ~ 10 10 cm -3 , and the average den- 
sity of the ionized region is lower than this. Furthermore, 
the ionization front is preceded by a photodissociation front 
whose extent always exceeds that of the ionized region, fur- 
ther ensuring that our ray-tracing scheme can safely ignore 
the presence of molecular species. 

The emission rate of H I, He I and He 11 ionizing photons 
is given by 



N Sm = 



^dv 
hu 



(3) 



where h is Planck's constant, ass the Stefan-Boltzmann con- 
stant, and z/ m m the minimum frequency required for ioniza- 
tion of the relevant species (H I or He 11). For simplicity we 
do not distinguish between the H 11 and He 11 regions. We 
assume the sink emits a blackbody spectrum B v with an 
effective temperature T e ff, which depends upon the evolving 

stellar radius and lu minosity. 

As described in IGreif et all (I2009D . the integral on the 
right-hand side of the ionization front equation is discretized 
by the following sum: 



n e n+r dr • 



riejn+^ri An 



(4) 



where Ar^ is the radial extent of each bin i, and the sum 
ranges from the sink particle to the current position of the 
Ffront. The left hand side of the ionization front equation 
is similarly discretized by: 



2 dn 
n " n I!- At 



(•>) 



where the sum now extends from the I-front position at the 
previous timestep to to its position at the current timestep 
to + At. The above-described ray-tracing is performed sep- 
arately for the H 11 and He ill regions. 

This ray-tracing scheme utilizes the simplifying 'on - 
the-spot' approximation (e.g. lOsterbrock fc Ferlandl I2006T J ■ 
which we will argue in Section 3.3 is a sufficiently accurate 
assumption. Our scheme also does not separately account for 
ionization by photons produced from recombination of He 11 
to He 1. Instead, Equation 2 includes He 11 in the n+ term. 
For simplicity, in terms of I-front evolution, He I and H I are 
effectively treated as the same species, as are He 11 and H II. 
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Figure 1. Evolution of various properties of the growing protostar according to our analytical model. Solid blue lines represent the 
protostellar values used in our 'with-fccdback' simulation. Dashed lines represent the ZAMS stellar values for the current sink mass. 
Dotted lines represent the 'slow contraction' case in which the accretion rate initially evolves in the same fashion as in the 'with-feedback' 
simulation, but then holds steady at 10 — 3 Mq yr . (a): Protostellar luminosity. Red dash-dotted line is the estimate for Lkh of a 15 
Mq star of radius 10 Rq. (b): Effective temperature, (c): Protostellar radius, (d): Ionizing luminosity, Ni on . Dashed-triplc-dottcd line 
represents the accretion rate of neutral particles onto the sink, extrapolated from the powerlaw fit to the sink accretion rate over the 
first 500 years (M oc t^cc 6 ' see section 3.5). Note how setting i?* = -RzAMS after 1000 yr in our simulation yields a good approximation 
to the more physically realistic 'slow-contraction' case. Both also predict a break-out of ionizing radiation beyond the sink at ~ 1000 yr, 
when iVi on exceeds the influx of neutral particles. 



Nevertheless, the helium abundance is small, ~ 0.08, and 
the typical He n abundance in the ionized region is slightly 
less than this value. Ionizing photons from He II recombi- 
nation will thus be < 0.1 of those from H n recombination, 
so the role of He II in the I-front evolution is relatively in- 
significant. For our star, which typically has T e g = 40,000 K 
(see Section 2.7), the He n and H II r egions indeed have the 
same radial extent as expected (e.g. lOsterbrock fc Ferlandl 
2006), while T c g is too low for an He III region to form. 

We note that we do not resolve higher-density regions 
within the sink, and thus do not directly simulate the disk 
self-shieldi ng and absorption of ionizing photons on this 
scale (e.g. iMcKee fc Tar]|2008r i. In order to model the ef- 
fect of self-shielding on sub-sink scales, we set the I-front 
radius to zero along all rays that encounter bins with den- 
sity greater than 5 x 10 9 cm -3 . This is the typical density 
of gas within the disk (see Section 3.1), and our prescription 



thus allows the directions along the dense molecular disk to 
be shielded while the more diffuse polar regions are the first 
to become ionized. This simple modeling of self-shielding 
also assumes that the large-scale disk surrounding the sink 
has the same orientation as the unresolved sub-sink disk. 

2.6 Photoionization and Heating 



Particles determined to be within the extent of the H 11 and 
He ill regions are given additional ionization and heating 
rates in the chemistry solver: 
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where F v is the incoming specific flux and cr„ the photo 
ionization cross section. For a blackbody, we have 



F, 



(8) 



4asBT e 4 ff r 2 

where r is the distance from the sink. 

Finally, H2 dissociation by LW radiation (11.2 to 13.6 
eV) is described by 

fe H2 = 1.1 x 10 8 / shicld F LW s" 1 (9) 

jAbel et al.lll997ft , where Flw denotes the radiation flux, in 
units of erg s" 1 cm -2 Hz -1 , at hV — 12.87eV, and /shield is 
the factor by which H2 self-shielding reduces the LW disso- 
ciation rate. This self-shielding factor depends upon the H2 
column density Nn 2 ■ With the above ray-tracing scheme, we 
determine Nu 2 along each ray by summi ng the contribution 
from e ach bin. We then use results from iDraine fc Bertoldil 
( 1996) to determine the value for /shi eld- We note a recent 
upda t e to their /shield fitting formula (I Wolcott- Green et al.l 
l201ll ; IWolcott-Green fc Haimanl l201ll ). but do not expect 
this to significantly affect the results for our particular case. 
Because of the unusually high column densities within the 
molecular disk (jV h, > 10 26 cm" 2 ), w e calcu late /shield using 
equation 37 from IDraine fc Bertoldil l| 1996ft , which is more 
accurate for large JVh 2 than their simple power-law expres- 
sion in their equation 36. 

These heating, ionization, and dissociation rates are ac- 
counted for at every timestep. They are updated every 10 
yr as they evolve with the protostellar mass and accretion 
rate, because these determine the values of L* and T c g. The 
ray-tracing procedure is perform e d eve ry fifth timestep. As 
discussed in IWhalen fc Normanl |2006;), it is important to 
update the ray-tracer often enough to correctly model the 
propagation speed of the I-front. We discuss a test later in 
Section 3.3 to show that our ray-tracing procedure is per- 
formed with sufficient frequency to correctly model the I- 
front growth. 

Our numerical method also neglects heating and ion- 
ization of gas ahead of the I-front by hard UV photons. 
Thus, instead of a gradual increase in temperature, ion- 
ized gas is instantly heate d to > 20, 000 K. As discussed in 
IWhalen fc Normanl <|200Sft and (2008b), this may artificially 
suppress I-front instabilities. This includes thin-shell insta- 
bilities driven by H2 cooling. These may be somewhat mit- 
igated by LW radiation, however, particularly at the early 
times studied in our simulation when the I-front remains 
in relatively close proximity to the star (< 0.1 pc) and H2 
shielding at the I-front edge is not too severe {Nu 2 ~ 10 16 
cm" 2 ). One may guess that such instabilities would allow 
ionizing radiation to be more easily channeled away from 
the stellar disk, allowing greater mass inflow onto the star 
(e.g. Whalen & Norman 2008b). Though currently too com- 
putationally expensive, a future study of 1-front instabilities 
on sub-parsec scales using multifrequency radiative transfer 
would be highly worthwhile. 



2.7 Radiation Pressure 

For simplicity, we do not include effects of radiation pressure 
in our calculation. We can estimate the effect of radiation 
pressure due to Thomson scattering by comparing the typ- 
ical protostellar luminosity (L* — 10 s Lq, Fig. [TJ with the 



Eddington luminosity: Z/Edd = 4:irGM,m p c/<jT, where M* 
is the stellar mass, m p is the mass of a proton, and <tt is 
the Thomson scattering cross section. For our typical stel- 
lar mass of 20 Mq , Z/Edd — 6 x 10 5 Lq . L* is several times 
smaller than this, so we would not expect Thomson scatter- 
ing to be dynamically important. 

Following simila r discussion in, e.g., lHaehneltl l|l995ft ; 
iMcKee fc Tanl (|2008ft , and given balance between ionizations 
and recombinations, we estimate the distance from the pro- 
tostar at which pressure from ionizing radiation will domi- 
nate over gravity by comparing radiative and effective grav- 
itational forces: 



hui 



= F, 



grav,cff 

> E dd GM,/im H 



(10) 



where 0Edd = 1 — Z/» /Z/Edd — 0.83. Using typical H 11 region 
densities of 10 7 cm -3 and hvi = 13.6 eV, Pion will become 
greater than Fg Ta ,v,ett beyond distances of ~ 100 AU, similar 
to the scale of the gravitational radius r g . However, beyond 
r g gas pressure should begin to dominate over both gravity 
and radiation pressure. Gas pressure is approximately given 
by nk B T > 3 x 10" 5 dyn cm" 2 within the 20,000 K H 11 
region. Ionization pressure is approximately 



N fhv. 



Aur 2 



(?) 



(11) 



which is ~ 3 x 10" 5 dyn cm" 2 for r = 200 AU and N = 



5 x 10 4 



Pi on will rapidly fall with radius, so on the 



scales which we study of 50 AU and beyond, gas pressure will 
usually dominate over ionizing radiation pressure. However, 
future work should include the effects of Pi on , particularly 
studies which examine very early I-front growth on scales 
below ~ 100 AU. 

Let us also estimate another potentially important ef- 
fect, radiation pressure from Lyman-a (Lya) photons. As- 
suming opaque conditions, pressure from Lya can be written 
aS Pet — — g U a = 47rJ a /3c where u a and J a are the energy 
density and mean intensity of the Lya radiation (see discus- 
sion injMcKe£_&Taii[|2008j) ; _A_^ is found in, 
e.g., lBithelll (|l990ft : IOh fc Haimanl l)2002ft . 



Pc, = 



Aivr 2 \ c 



TLya 



(12) 



where Nhya is the rate of emission of Lya photons by the 
star and s T is the Lya photon path length in units of the op- 
tical depth of the region of gas in quest ion, tlvq (see also dis - 
cussion and more de tailed equation in lMcKee fc Tan1l2008ft . 
From equation 21 of iBithelll l|l990ft . we estimate s r /rL yQ1 to 
be ~ 5000 for the disk gas and ~ 100 for the ionized region. 
However , this may be further reduced by the motion of th e 
gas (e.g. lBithelllll99d : lHaehneitlll995l : lMcKee fc Tanll2"008h . 
The Lya photons thus travel much more freely through 



the ionized region. Estimating A^ya 



N 



5 x 10 4 



for the ~ 20 Mq protostar, we then find upper limits of 
P„-5x 10" 3 dyn cm~ 2 at the edge of the 1000 AU disk, 
and 10~ 5 dyn cm" 2 at the e dge of the 10 4 AU H 11 region. 

As discussed in detail in lMcKee fc Tar] l|2008ft . Lya ra- 
diation pressure can be estimated to significantly slow pro- 
tostellar accretion when it is over twice the ram pressure of 
infalling gas, P a > 2pvg, where vg is the free-fall velocity of 
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infalling gas. For gas in a rotating disk, the requirement be- 
comes P a > p^Kcp — At the disk edge where n ~ 10 9 



10~ 3 dyn cm -2 . At the edge of the H n 
region where n ~ 10 7 cm -3 , we have 1pv% < 10~ 5 dyn 
cm~ 2 . These values are indeed exceeded by the correspond- 
ing strength of P a . In contrast, within a few hundred AU 
of the main sink where n ~ 10 10 - 10 11 cm 3 , P a is up to 
several times smaller than pi>K C p, and the gravitational force 

of the sink is more dominant. 

Consistent with the findings of iMcKee fc Tanl (|2008h . 
Lya may indeed break out when a star has reached 20 Mq, 
or upon reaching even larger masses for the case of gas 
undergoing less rotation. Lya pressure is thus expected to 
break out only after the star has become massive enough to 
develop an H II region. We would furthermore expect Lya to 
break out in the polar regions before affecting the disk plane, 
as the polar regions are more diffuse. Afterwards, Lya will 
easily escape through the polar cavities, relieving pressure 
closer to the disk plane through which the protostars ac- 
crete. We thus do not expect that inclusion of Lya pressure 
would significantly change our results, though future work 
should confirm this through more detailed simulations. 



2.8 Protostellar evolution model 

Our ray-tracing algorithm requires an input of protostellar 
effective temperature T e s and luminosity L». This is the sum 
of Lace, the accretion luminosity, and Lphoto, the luminosity 
generated at the protostellar surface: 



+ L 



photo 



GM,M 



+ L 



photo i 



(13) 



where M, is the protostellar mass, M the accretion rate, 
and R* the protostellar radius. The photospheric luminos- 
ity is either due to Kelvin Helmholtz (KH) contraction, or 
due to hydrogen burning which begins once the protostar 
has reached the zero-age main sequence (ZAMS). For sim- 
plicity, we set Lphoto — Lzams, thus assuming a robust 
upper limit in the luminosity at later phases in the pro- 
tostellar evolution. Initially, however, we assume no contri- 
bution from Lphoto until -R* reaches the ZAMS radius. This 
is a reasonable assumption, pa rticularly given the valu es of 
Lphoto as determined by, e.g., iHosokawa et al.l (|2010h . We 
find that L acc should be much greater or similar in magni- 
tude to Lphoto until the ZAMS is reached in our model (see 
Fig. [TJ . Our assumption is also robust in view of uncertain- 
ties regarding h o w ra pidly KH contraction proceeds. While 
IHosokawa et alj (|201fj| ) find that a massive (~ 10 Mq) ac- 
creting protostar will typically contract to the ZAMS radius 
on timescales of ~ 10 4 yr, in our model we set our protostar 
to the ZAMS values much earlier. However, the luminos- 
ity during KH contraction is indeed similar in magnitude t o 
Lzams for a star of the same mass (|Hosokawa et al.ll2010D . 
Nevertheless, it is important that future simulations self- 
consistently couple protostellar evolution and accretion flow 

under fee dback. 

As in lSchaller et al.1 l|l992l) and IHosokawa et alj i|2010t ). 
we determine Lzams using a simple fit to the stellar mass: 



Lzams = 1.4 x 10 L G 



M* 

10Mr. 



(14) 



Each time L* is updated, we assume M* is equal to the 
sink mass. Due to the discrete nature of sink accretion in an 
SPH simulation, instead of calculating M at each timestep, 
we determine M by averaging the sink mass growth over the 
previous 100 yr, updating M every 10 yr. 

We es ti mate R, in the same fashion as in 
IStacv et all J2010fl. wh i ch wa s based on the prescrip- 
tion of lOmukai fc Pallal l |2003t ). We find that during the 
adiabatic accretion phase, -R* grows as 



R, 



50R e 



M, 
M© 



1/3 



M 



1/3 



(15) 



where Mud — 4.4 x lO^Mgyr -1 is a fiducial rate, typical 
for Pop III accretion. Throughout this phase, we assume 
there is not yet any contribution from Lphoto- During the 
subsequent phase of KH contraction, the radius will shrink 
according to 



R* 



MORp 



/ M \ ( M, 



\MudJ V 10M e 



(16) 



We estimate that the transition from adiabatic accretion to 
KH contraction occurs when the value of R,n falls below 
that of R,i. During this phase, our model again assumes 
no luminosity contribution from Lphoto, and that L acc is the 
main contribution to the luminosity. KH contraction will 
halt once the star has reached the ZAMS, at which point we 
set R* equal to the ZAMS radius, 



-Rzams = 3.9R(] 



AL 



(17) 



y 10M Q 

(e.g. IHosokawa et al.ll2010l ). We set R* equal to -Rzams when 
the value for R,n falls below -Rzams- 

If the calculated accretion rate drops to near zero, 
then the radial values for the adiabatic and KH contraction 
phases will become vanishingly small. If this occurs before 
the sink has been accreting for a KH time and reached the 
ZAMS, the accretion rate is set to the previous non-zero 
value in order to get more realistic values for R*i and R*n. 
This allows us to avoid setting R, = -Rzams too early in the 
protostar's evolution. If, however, the accretion slows after 
the sink has been in place for more than its KH time, we 
assume the star has reached its ZAMS radius, and we set 
Lacc — 0, R t — -Rzams, and L* = Lzams- Note that typical 
KH times, where tKH = GM% / R*L*, range from 1000 yr 
for a large and rapidl y accreting 10 Mq protostar (see e.g. 
IHosokawa et~ai1l2010l ) to ~ 4 x 10 4 yr for a 15 M mam 
sequence star. The typical KH luminosity for a 15 Mq star 
is L KH ~ 5 x 10 4 L© (see Fig.[TJ|. 

Given our averaging scheme in which a minimum of 
one 0.015 M gas particle can be accreted over 100 years, 
this gives an effective minimum measurable accretion rate 
of 1.5 x 10~ 4 M Q yr" 1 . However, for M* > 10 M Q , this 
minimum accretion rate still yields a value of that is 

smaller than -Rzams- In this case, we again set the proto- 
stellar luminosity and radius to its ZAMS values once the 
accretion rate has dropped to < 10 -4 Mq yr -1 . In our case, 
the measured accretion rate drops very quickly after 500 
yr. At this point the star has reached 15 Mq, is still un- 
dergoing adiabatic expansion, and has txH ~ 1000 yr. The 
star then begins rapid KH contraction until the measured 
accretion rate drops below 10 -4 Mq yr -1 at > 1000 yr, 
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-2000 2000 4000 -2000 2000 4000 

time [yr] time [yr] 



Figure 2. Evolution of various disk properties. Time is measured from the point at which the first sink particle forms. In each case, solid 
lines represent the 'no-feedback' case while the dotted lines denote the 'with-feedback' case, (a): The average density of the gas within the 
disk, (b): The ratio of the average radial to the average rotational velocity of the disk particles, Xrad- ( c ) : The average temperature of the 
disk particles, (d): The average angular velocity of the disk gas with respect to the disk center of mass, Q. Note how the 'with-feedback' 
and 'no-feedback' cases diverge after 1000 yr. 



not rising above this value again until ~ 2000 yr (see in- 
crease in luminosity at this time in Fig. [TJ . Though within 
the simulation we set R* = -Rzams as soon as the averaged 
accretion rate falls below 10~ 4 Mq yr -1 , in reality the pro- 
tostar is better described by a more gradual approach to 
-Rzams • In Figure [TJ we show the protostellar values used in 
the simulation along with a more realistic 'slow-contraction' 
model which follows the same accretion history as the 'with- 
feedback' case until reaching an asymptotic growth rate of 
10~ 3 Mq yr -1 . The 'slow-contraction' model is then held 
at this rate, which is similar to the fiducial value used in 
Equation 11, and is also the typical accretion rate found at 
late times in our 'no-feedback c a se' as well a s the s imula- 
tions of, e.g., iGreif et all (|201lh : ISmith et all (|201lh . This 
model well-matches the prescription used in the simulation, 
particularly in effective temperature and ionizing luminos- 
ity, and both predict that ionizing radiation will exceed the 
influx of neutral particles, and that break-out beyond the 
sink occurs at ~ 1000 yr. Though the protostar most likely 
does not reach the ZAMS within the time of our simulation, 



our simple model serves as a reasonable approximation for 
any unresolved accretion luminosity. 

We also note that iHosokawa et al.l (|2010l ) find that, for 
a given accretion rate, primordial protostars undergoing disk 
accretion will have considerably smaller radii than those ac- 
creting mass in a spherical geometry. The rapid contraction 
of the protostar between 500 and 1000 years therefore serves 
as an idealized representation of the sub-sink material evolv- 
ing from a spherical to a disk geometry as it gains angular 
momentum. We emphasize, however, that the unrealistically 
rapid contraction to the ZAMS in our model, along with the 
incomplete inner-disk self-shielding, leads to an overestimate 
in the strength of feedback. In this way we underestimate 
the mass of the star. In addition, the recent work by Smith 
et al. (2011b) found that Pop III protostars undergoing vari- 
able accretion rates may have very large radii (100-200 Rq) 
for most of their pre-main sequence lifetimes. This again in- 
dicates a probable underestimate of the protostellar radius 
and an overestimate of T c g in our model. However, note 
that Smith et al. (2011b) assumed a ' hot' spherical accre- 
tion model (e.g. IHosokawa et al.l [2010. 2011b), which will 
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time [yr] time [yr] 

Figure 3. Evolution of disk mass over time. Solid lines are for the 'no-feedback' case. Dotted lines represent the 'with-feedback' case. 
(a): Black lines show the mass of the total star-disk system. Disk mass (green lines) is taken as dense (n > 10 9 cm — 3 ) gas with an 
H2 fraction greater than 10 -3 . (b): Red lines show the mass of hotter, non-molecular material in the same density range as the disk, 
and blue lines denote the total mass of the outer envelope, which is defined as comprising all sinks and n > 10 s cm — 3 gas. Note how 
radiative feedback causes the total disk and envelope mass to decline after approximately 2000 yr. N-body dynamics in the 'no-feedback' 
case causes an even steeper initial decline in disk mass compared to radiation in the 'with-feedback' case. 



yield significantly larger radii than the 'cold' thin disk ac- 
cretion model. Though disk thickness as well as high disk 
temperatures and accretion rates may render 'hot' accretion 
the more realistic case, the true radial evolution of the star 
is likely to be in between these two possibilities. 



As discussed in ISmith et al.l (|201l[ ). the sink particle 
method requires several simplifying assumptions when con- 
structing the protostellar model. By setting the protostel- 
lar mass equal to the sink mass, we are assuming that the 
small-scale disk which likely exists within the sink region 
has low mass compared to the protostar. We also assume 
that the accretion rate at the sink edge is the same as the 
accretion rate onto the star, when in reality after gas enters 
within r acc it must likely be processed through a small-scale 
disk before being incorporated onto the star. However, our 
assumption may still be a good approximation of physical 
reality, given that the primordial protostellar disk study of 
Clark et al. (2011b) implies that the thin disk approx imation 
woul d probably not be valid on sub-sink scales (e.g. IPringld 
Il98ll ). and that strong gravitational torques can quickly 
dr ive mass ont o the s tar. Furthermore, as also pointed out 
in ISmith et al.l (|201ll ). our averaging of the accretion rate 
over a number of timesteps serves to mimic the buffering of 
accreted material by the sub-sink disk. The inputs to our 
protostellar model are thus necessarily approximate, as is 
the protostellar model itself. Nevertheless, it is sufficient to 
provide an exploratory picture of how ionizing feedback will 
affect Pop III accretion within a cosmological setup. 



3 RESULTS 

3.1 Disk Evolution 

The gas in both the 'no-feedback' and 'with-feedback' cases 
underwent disk formation, fragmentation, and the emer- 
gence of several sinks. We show the evolution of various disk 
properties in Figure [2] Because of the imprecision involved 
in determining which gas particles comprise the disk, for 
simplicity we define the disk as consisting of particles with 
number density greater than 10 9 cm -3 and with an H2 frac- 
tion greater than 10 -3 . This way the disk only contains cool 
and dense gas that has not been subject to ionization or 
significant H2 destruction. From Figure [3] we see that the 
disk structure is growing in mass well before the first sink 
forms. This central gas is already rotationally dominated, as 
indicated in Figure [2] by the low values of Xrad, which is the 
average radial velocity of the gas particles divided by their 
average rotational velocity, Xrad = Urad/frot- Velocities are 
measured with respect to the center of mass of the disk. 



3.1.1 No-feedback case 

After sink formation, the 'no-feedback' disk steadily grows 
in mass (Fig. [3j as angular momentum causes it to expand 
in radius and become somewhat lower in average density 
(panel a in Fig. [5} . The disk growth is halted at nearly 2000 
yr due to its gradual disruption through N-body dynamics 
of the sinks. One of the sinks is ejected at ~ 500 yr, after 
growing to only ~ 1 M©. The ejection occurs immediately 
following the merger of the two other sinks, at a time when 
the three sinks are close together in the center of the disk 
and subject to N-body dynamics. The sink accretes no more 
mass after its ejection. It initially moves in a direction per- 
pendicular to the disk plane at ~ 5 km s" 1 with respect 
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to the disk center of mass, and the maximum distance be- 
tween the two sinks approaches 3000 AU at approximately 
2500 yr. This increase in distance between the two sinks is 
mostly due to the motion of the main sink, however. In this 
N-body interaction, the main sink gains a larger velocity, 
initially moving with respect to the disk center of mass at 
~ 10 km s _1 . It travels parallel to the disk plane and pulls 
the disk along with it. The rapid motion of the main sink 
disrupts the high-density gas, which transforms from a disk 
structure to a more diffuse tidal tail, eventually causing the 
total measured disk mass to slightly decrease (Fig. [3]). The 
rotational structure is also disturbed, as indicated by a peak 
in Xrad when the sink is ejected. After approximately 3000 
yr, the tidal tail begins to recompress, causing the dense gas 
to be more dominated by radially inward motion. This is in- 
dicated by the increase of Xrad to ~ 1 as well as an increase 
in the average disk density (Fig. 

From the solid red line in Figure [3J we also note a 
smooth early growth of hot dense gas over the first 1000 
yr. This begins almost as soon as the first sink is formed, 
and the main sink thus provides an early source of heating. 
This occurs as the gravitational potential of the sink, which 
grows to > 10 Mq in only 200 yr, heats the surrounding gas 
up to ~ 7000 K, or c s ~ 7 km s~ (Fig. 2}, corresponding 
to approximately the virial temperature of the sink: 

rt ,. Gto ^ 1()4K (lg) 



The mass of hot dense gas undergoes an equally smooth 
decline over the next 1000 yr in the 'no-feedback' case. A 
minimal amount of dense gas is newly heated after 1000 yr 
because there is no longer a dense disk structure entirely 
surrounding the main sink. Most of the dense gas instead 
trails behind the sink in the tidal tail, and the gravitational 
heating provided by the sink is now imparted directly to 
the lower-density particles. Sink gravitational heating thus 
halts the growth of the outer envelope, which we define to 
include all gas particles at densities above 10 8 cm -3 as well 
as the sinks (solid blue line in Fig. [3]). The 'heat wave' travels 
beyond the disk at approximately the sound speed, causing 
the decline in envelope mass beginning at around 2500 yr in 
Figure [3] and reaching distances of > 7000 AU (see Figs. 2] 
and[5]) by the end of the simulation. At this time, heated gas 
thus resides almost entirely at low densities (n < 10 9 cm -3 ) 

This d i sk an d envelope evolution resembles that in 
IStacv et al.l (|2010h . particularly in that early gravitational 
heating by the sink caused the developme nt of a phase of hot 
and dense gas. However, the hot phase in lStacv et al.l (|2010l ) 
fell mostly between densities of 10 8 and 10 12 cm -3 , and the 
mass of both the envelope and total star-disk system showed 
a steady increase throughout the simulation. In contrast, in 
our current 'no-feedback' case, the hot phase lies mostly 
between densities of 10 6 and 10 9 cm -3 , and both the disk 
and envelope decline in mass after 2000 yr. This difference 
between simulations ultimately arises from the statistical 
variation in sink N-body dynamics, and the corresponding 
response from the surrounding gas. 



3.1.2 With-feedback case 

The 'with-feedback' case also exhibits a peak in Xrad coin- 
cident with the period of rapid sink formation, with a max- 



imum of three sinks in the disk at any given time. This 
causes the disk to be dominated by N-body dynamics and 
disrupts the rotational structure, but without ejection of any 
sink. However, the disk growth is soon slowed by a different 
mechanism than that in the 'no-feedback' case. Just as in 
the latter, we also note a smooth early growth in the mass of 
hot dense gas (Fig.[3j, initially sourced by the gravitational 
potential of the main sink. Once the sink is large enough to 
emit significant amounts of LW and ionizing radiation, after 
~ 900 yr of accretion, infalling mass that would otherwise 
be incorporated into the disk is instead heated to become 
a hot neutral shell of several thousand Kelvin enclosing a 
smaller ionized bubble around the disk (see bottom panels 
of Fig. [3] and Fig. [5]). This leads to the continued increase 
in the mass of the hot dense gas apparent in Figures [3J and 
[6] beyond that seen in the 'no-feedback' case. As the ion- 
ization front and hot neutral region expand and rarefy the 
gas, both the disk and hot dense regions lose mass. This is 
visible as a drop in mass of the hot dense gas at 1500 yr 
in Figure [3J once the pressure wave has reached beyond the 
dense region into n < 10 9 cm -3 gas. Throughout this early 
disk evolution, the total mass of the outer envelope steadily 
increases (blue line in Fig. [3J until the same pressure wave 
passes through the edge of the envelope after 2000 yr, after 
which both the envelope and disk gradually lose mass until 
the end of the simulation. 

At the same time, angular momentum causes the disk 
to expand in radius, making the disk more diffuse (see panel 
a of Fig. |5J. Once the hot pressure wave travels into the 
outer envelope, the remaining disk mass is able to level off 
in density after 2000 yr. Meanwhile, the shielded regions of 
the disk remain steady in their rotational structure (panels 
b and d of Figure [2]). It is interesting to note that, because 
of disk shielding, for the first 4000 yr N-body dynamics had 
a greater effect on the disk growth in the 'no-feedback' case 
than protostellar radiation in the 'with-feedback' case. How- 
ever, the scale of radiative feedback is much larger than that 
of N-body dynamics, as can be seen by the greater decline 
in mass of the outer envelope in the 'with-feedback' case. 



3.2 Fragmentation 

Consistent with ISmith et al.l |201ll ). the accretion luminos- 
ity does not heat the disk sufficiently to prevent fragmen- 
tation. Radiative feedback does slightly lower the overall 
fragmentation rate, however, as a total of eight sinks are 
formed in the 'no-feedback' case versus five sinks in the 
'with-feedback' case. In both cases the second sink forms 
only ~ 100 yr after the initial sink, but, similar to most of the 
sinks formed, it is quickly lost to a merger. The next sinks 
to form and survive to the end of the simulation are created 
at 300 and 200 yr in the 'no-feedback' and 'with-feedback' 
cases, respectively. This very quick fragmentatio n is similar 
to wh at wa s described in Clar k et al (2011b), iGreif et al.l 
( |201lt ). and lSmith et all l|201lf ). The last sink forms as late 
as 2000 yr in the 'no-feedback' simulation (see Table [!}. 

To understand the formation of the initial disk gravi- 
tational instability, we first check to see if the disk satisfies 
the Toomre criterion: 



ttGE 



< 1 , 



(19) 
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Figure 4. Temperature versus number density for both cases at various times in the simulations, (a): 'No-feedback' case at 2500 yr. (b): 
'No-feedback' case at 5000 yr. (c): 'With-feedback' case at 2000 yr. (d): 'With-feedback' case at 3000 yr. Note how, in the 'no-feedback' 
case, there is only a light stream of particles with n > 10 9 cm -3 that is accreting onto the main sink. The gravitational potential well 
of the main sink leads to the heating of a growing region of lower-density gas. In the 'with-feedback' case, there is an expanding ionized 
region with temperature of 20,000 K along with a larger region of hot neutral gas with temperature of several thousand Kelvin. 



where c s is the soundspeed, E the disk surface density, and 
k the epicyclic frequency, which for Keplerian rotation is 
equal to the disk angular velocity Q. For the first few 100 
yr both disks have average temperatures of around 1000 
K, so c s ~ 2 km s _1 . A/disk — 25 Mq, and the disk ra- 
dius is nearly 1000 AU, yielding a disk surface density of 
E ~ Mdisk/i"-Rdisk ~ 70 g cm -2 . We also approximate k to 
be 3x 10 -11 s _1 , so Q ~ 0.4, satisfying the Toomre criterion. 

The ability of the disk to reach Toomre instability is 
aided by the rapid cooling time of it s gas. The criterion for 
fragmentation described in iGammiel (|200ll ) is i coo i < 3fi _1 , 
where t CO oi is the cooling time of the disk gas. We find that 
the typical value of r coo i is ~ 50 yr, while 3f2~ ~ 3000 yr, 
easily satisfyin g this criterion, also refe r red to as the 'Viscous 
Criterion' bv lKratter fc Murrav-Clavl (|201ll ). 

Figures [S] and [5] show the density and temperature mor- 
phology within the central 10,000 AU. The multiple sinks 
and clumpy disk structure are easily visible here, as expected 
from the low Toomre parameter, though the particular shape 
of the disk in each case is very different. The 'no-feedback' 
case has a clearly visible bifurcated temperature structure, 



sink 


*form [yr] 


A/final [M ] 


r-i„it [AU] 


r&nal [AU] 


1 





27 








2 


300 


0.9 


100 


2330 


3 


2000 


2.75 


70 


83 



Table 1. Formation times, final masses, distances from the main 
sink upon initial formation, and distances from the main sink at 
the final simulation output in the 'no-feedback' case. We include 
the sinks still present at the end of the simulation (5000 yr). 

with a cool disk and surrounded by warm gas. The 'with- 
feedback' also has cool disk gas, but the central region is 
much more dominated by an hour-glass shaped bubble of 

hot gas. 

As discussed in iKratter et all <|2010h . the actual mini- 
mum value of Q can vary with disk properties such as scale 
height, and from Figures [5] and [6] we see this varies through- 
out the disk evolution for both cases. Thus, it is also nec- 
essary to consider other disk properties, particularly the in- 
fall rate of mass onto the disk. The numerical experiments 
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Figure 5. Projected density and temperature structure of central 10,000 AU without protostellar feedback at 1000 yr (left column), 
2000 yr (middle column), and 5000 yr (right column) after initial sink formation. Asterisks denote the location of the most massive sink. 
Crosses show the location of the ejected sink. Diamonds mark the locations of the other sinks. Top row : Density structure of the central 
region in the x-y plane. Second row : Density structure of the central region in the x-z plane. Third row : Temperature structure of the 
central region in the x-y plane. Bottom row : Temperature structure of the central region in the x-z plane. Note the formation of a cool 
disk and tidal tail structure along with the growth of a surrounding bubble of warm gas. The ejected sink (cross) has a mass of 1 Mq. 
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Figure 6. Projected density and temperature structure of gas under LW and ionization feedback at 1000, 2000, and 3000 yr after initial 
sink formation. Asterisks denote the location of the most massive sink. Crosses show the location of the second-most massive sink. 
Diamonds are the locations of the other sinks. For the rows and columns, we adopt the same convention as employed in Fig. [5] Note the 
growth of a roughly hour-glass shaped structure of hot gas surrounding the disk as the I- front expands into the low-density regions. The 
hot region expands well beyond the disk by 2000 yr. 
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sink 


iform [yr] 


M flna i [M Q ] r init [AU] 


r fiaa l [AU] 


1 





19 





2 


200 


9.4 110 


440 



Table 2. Same as Table 1, but for sinks remaining at the end of 
the 'with-feedback' case at 4500 yr. 



of iKratter efall (|2010h determined that further fragmen- 
tation will occur if the mass infall rate onto the disk is 
sufficiently high, M- lTL > c^/G, such that the disk can no 
longer process the new material quickly enough. The star- 
disk system in the 'no-feedback' case grows at approximately 
5 x lO" 3 M yr- 1 over the first 2000 yr of infall (Fig.©. For 
the 'with-feedback' case the star-disk system grows at this 
same rate, but for a few hundred years longer. Comparing 
this to the average gas soundspeed in the disks, ~ 2 km s _1 
in both cases, we find that Mi n is almost three times greater 
than c 3 /G. Thus, as expected, the disks in both cases frag- 
ment during their initial growth phases, though this phase is 
slightly more unstable for the cooler gas of the 'no-feedback' 
case (panel c of Figure(2]|, for which the formation of the last 
sink also coincides with a final 'bump' in the disk mass and 
density at 2000 yr. Though there is a similar bump in disk 
mass at this time for the 'with-feedback' case, the gas is a 
few hundred Kelvin warmer. This, along with the steadily 
decreasing average density of the disk (panel a of Figure [2]) , 
is sufficient to keep the disk stable at later times. No new 
sinks form once the growth of the disks is halted. 

Once the disk has become unstable to fragmentation, 
the continued gravitational collapse of the disk fragments 
furthermore requires that the gas cool quickly enough to 
overcome opposing effects such as pressure and tidal forces 
(t he 'Stalling Criterion'), a s des cribed in the recent work 
of IKratter fc Murrav-Clavl (|201lf ). They find that this cri- 
terion should be easily met for molecular gas (7 = 7/5), 
which applies to the densest regions of the disks in our sim- 
ulations. As mentioned above, for the first few hundred yr 
after sink formation, the average cooling time of the disk is 
fc ooi ~ 50 yr. We compare this to the /3 factor described 
in IKratter fc Murrav-Clavl (|201lh to find /3 = r coo iri < 0.1. 
The critical value of /3 necessary for free-fall collapse of the 
fragment ranges from 2.5 to 13 depending on 7, and the 
disks easily satisfy this criterion, thus leading to the early 
formation of sinks within the disks. 

However, in the densest regions of the disk, the 



' Collis ional Criterion' described in IKratter fc Murrav-Clavl 
(2011), which requires that the sinks collapse on roughly 
the orbital timescale, is harder for these sinks to satisfy. 
When the disks are undergoing fragmentation during the 
first few hundred yr, the sinks form within ~ 100 AU of 
each other and orbit at approximately 5-10 km s _1 , yield- 
ing an orbital timescale on the order of 100 yr. This is not 
significantly longer than the free-fall timescale of the sinks, 
which is < 100 yr. This leads to sink merging before an en- 
tire orbit is complete, particularly during early times in the 
'no-feedback' case. The one merger that occurs for the main 
sink in the 'with-feedback' case occurs after multiple orbits, 
however, and is a result of migration. 



3.3 Evolution of Ionization Front 

The I-front initially appears around the main sink ~ 1000 
yr after the sink first forms, once it has grown to ~ 15 Mq 
(see panel d of Fig. [TJ. The morphology of the growing I- 
front is shown in Figure [7] where we can see that the 1- front 
expands as an hour-glass shape above and below the disk. 
This same morphology describes the growing neutral region 
as well (Fig. [8}, which we will discuss in the next section 
(3.4). As the I-front expands and widens to encompass a 
larger angular region, it dissipates the more diffuse disk gas 
above and below the mid-plane, leading to a decline in the 
disk scale height. Figure [5] compares the I-front evolution 
with tha t predicted by t he analy tical Shu champagn e flow 
solution |Shu et al.ll2002l . see also lAlvarez et al.l l2006). Dif- 
ferent analytical solutions can be found for the evolution 
of gas under the propagation of an I-front into a powerlaw 
density profile, which in our case is approximately p oc r~ 2 . 
The ratio of the un-ionized gas temperature to that of the 
H 11 region must also be specified for the analytical solution, 
and in our case we choose a ratio of 20,000 K to 1000 K, or 
0.05. The propagation of the I-front can then be described, 
assuming D-type conditions and that the I-front closely fol- 
lows the preceding shock, by a velocity of 



while the size of the I-front is 



(20) 



(21) 



where c s is the ionized gas soundspeed and x s is the position 
of the shock in similarity coordinates, which for our case is 
x s = 2.54 (see panel b of Figure |9}. 

Note that the typical neutral hydrogen abundance /hi is 
around 10 -2 in the ionized region, while typical densities in 
the latter part of the simulation are 10 7 cm -3 (panels c and 
d of Fig. |9]). For a hydrogen photoionization cross section of 
fion = 6 x 10 -18 cm 2 , this yields a mean-free-path of 

2 x 10 17 cm 



^mfp — 



1 



o-ionn m / H m[cm 3 ] 



~ 2 x 10 13 cm 



(22) 



where nm is the number density of neutral hydrogen. The 
value of Z m fp will usually be ~ 1 AU. This is much smaller 
than the typical length of a radial bin in the simulation, 
with the exception of the small 1.5 AU bins within 75 AU of 
the star, and also much smaller than the resolution length of 
the simulation. Recombination photons should therefore be 
reabsorbed before exiting their ray-tracing bin and entering 
neighboring bins. Our ray-tracing scheme thus loses little 
accuracy by applying the 'on-the-spot' approximation and 
ignoring the effects of diffuse radiation within the H 11 region. 

Multiple factors cause deviations from the analytical 
solution. For instance, the density structure of the gas is 
not spherically symmetric. Furthermore, the initially close 
proximity of the ionization front to the sink causes the grav- 
ity of the sink to have a non-negligible effect on the early 
H 11 region dynamics, a factor that is neglected in the an- 
alytical solution. However, this effect loses importance as 
the H 11 region grows well beyond the gravitational radius, 
r g ~ GM,/c s ~ 50 AU for a 15M Q star and 20,000K H 11 
region. The analytical solution also does not take into ac- 
count the continued infall of gas onto the central region. 
The outer envelope beyond the disk, which we take as the 
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Figure 7. Projected ionization structure of gas at 1500, 2000, and 3000 yr after initial sink formation. White lines depict the density 
contours of the disk at densities of 10 7 5 , 10 s , 10 85 , and 10 9 cm -3 . Box length is 40,000 AU. Note the pronounced hour-glass morphology 
of the developing ultra-compact H II region, roughly perpendicular to the disk. This structure gradually expands and dissipates the disk 
gas from above and below, reducing the scale height of the disk. 
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Figure 8. Projected H2 fraction of gas at 1500, 2000, and 3000 yr after initial sink formation. Gray lines depict the density contours 
of the disk at densities of 10 75 , 10 s , 10 8 ' 5 , and 10 9 cm . Box length is 40,000 AU. Note how the morphology of the neutral region is 
similar to that of the ionized region. The molecular inflow gradually becomes confined to the disk plane. 



gas with density greater than 10 8 cm 3 , grows for the first 
2000 yr at a rate of ~ 10~ 2 M yr" 1 , or 5 x 10 47 neutral 
particles per second. In comparison, the sink typically emits 
> 10 48 photons per second. Though the infall is not large 
enough to quench the H n region entirely, it is a large enough 
fraction of the ionization rate to cause the total mass of the 
H 11 region to fluctuate. The motion of the sink through 
the disk causes the immediately surrounding density struc- 
ture to vary as well, also contributing to the H 11 region's 
unsteady growth. Finally, we note that the peak in the I- 
front radius at ~ 2500 yr (panel b of Figure [§]) corresponds 
to an increase in the luminosity of the main sink due to a 
concurrent enhancement of the accretion rate. 

This fluctuation is similar to what lGalvan-Madrid et ail 

(2011) described in their numerical study of hypercompact 
H 11 regions. However, our Ffront evo lution differs signif- 
icantly from the analytical study by lOmukai Sz Inutsukal 
(2002), where they found that a spherical free- falling en- 
velope would not be unbound by an I-front typically until 
the Pop III star reaches well over 100 Mq . In our study the 
gas does become unbound in the regions polar to the disk, 
thus indicating how variations in three-dimensional struc- 



ture and accretion flow play a crucial role in a Pop III star's 
accretion history. 

Because we do not resolve the gas on scales smaller than 
~ 50 AU, there may be unresolved substructure which would 
have provided shielding and altered the H 11 region growth. 
Our simulation only approximates this shielding effect with 
a simple prescription (Section 2.5) that assumes a smooth 
and dense sub-sink disk that is coplanar with the large-scale 
disk resolved in our simulation. The true shielding within the 
sink may vary from this assumption, however, depending on 
the details of the sub-sink str ucture. 

Clark et al. (2011b) and iGreif et alj (|201ll ) find that 
disk formation indeed continues down to < 50 AU scales, 
though the disk fragments and clumps on these scales in- 
stead of maintaining a smooth structure. Were sub-sink 
scales in our simulation to contain a similar clumpy disk 
structure that was coplanar with our ~ 1000 AU disk, this 
would still shield much of the gas lying within the disk 
plane, as provided by our prescription. However, assuming 
no change in the available sub-sink disk mass, dumpiness 
in the disk will likely increase the escape fraction of ioniz- 
ing radiation from the sink, since extra radiation may leak 
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through 'holes' in the disk (see, e.g. I Wood fc Loebl feoOO). 
The growth of the H n region in our simulation would then 
be an underestimate. On the other hand, dumpiness along 
lines-of-sight perpendicular to the disk plane would provide 
extra shielding in these directions and slow the initial growth 
of the H II region in the polar directions. Our calculation 
would then have overestimated the H II region size, partic- 
ularly given the n 2 dependence of the recombination rate. 
Our shielding prescription is thus a median between these 
two possibilities. 

As for structure within the H II region, our prescrip- 
tion ensures that gas is not ionized unless its density is ap- 
proximately 100 times lower than our maximum resolvable 
density (~ 10 10 cm -3 versus 10 12 cm -3 ). Once ionized, gas 
particles may condense to 100 times their initial density be- 
fore structure within the H II region becomes unresolved. 
Formation of structure on small ~ 50 AU scales within the 
H II region does not occur in our calculation, however, and 
H II region substructure should not significantly alter its 
evolution. 

To test whether performing the ray tracing procedure 
every five timesteps was sufficient to accurately follow the 
growth of the I-front, we also performed a test simulation 
in which the ray tracer was updated at every timestep. The 
test simulation was initialized using the output of the 'with- 
feedback' simulation at 1000 yr, and then followed for the 
following ~ 1000 yr. As seen in panels a and b of Figure 
[9] increasing the frequency of the ray-tracing update led to 
little change in the I-front evolution, so updating every five 
timesteps was sufficient. 

3.4 Evolution of Hot Neutral Region 

Along with the H II region, suppression of H2 cooling by LW 
radiation, combined with continued gravitational heating by 
the sink, leads to the development of a growing region of hot 
neutral gas that expands in a bubble around the disk. There 
is a pressure wave of warm gas, sourced by the combination 
of the sink's gravitational potential and LW radiation, that 
expands at the soundspeed of approximately 7 km s _1 (see 
also Section 3.1). After 3000 yr this pressure wave has ex- 
panded to a radial extent of r pros ~ 5000 AU, matching well 
the predicted size of r prcs = c s t = 4400 AU. This corre- 
sponds to the region of gas in Figure [4] centered around 10 8 
cm -3 and sitting between approximately 3000 and 7000 K. 

This inner pressure wave is not the sole source of hot 
neutral gas, however. The hot neutral region also has an 
extended 'tail' of 7000 K gas that can be seen in Figure [4] 
This corresponds to the extended region that became ionized 
during the main sink's enhanced accretion phase, but then 
recombined and dropped to temperatures of 7000 K when 
the sink's luminosity declined again. LW radiation prevents 
this 'tail' of hot neutral gas from cooling back down again. 
Thus, LW radiation, gravitational heating by the sink, and 
recombination of ionized gas all contribute to the growth of 
the hot neutral region that encompasses the I-front. 

We note that the photodissociation front extends even 
further than this, beyond 10 5 AU (nearly a pc) from the sink 
by the end of the simulation. The edge of the photodissoci- 
ation front is comprised of neutral but cool gas. This pho- 
todissociation region can expand much more quickly than 
the gas soundspeed, and thus is much larger than the scale 



encom pass ed by the pr e ssure wave. Following I Johnson et al.l 
(|2007f ) and I Abel et all <|l997l ). we can estimate the dissoca- 
tion time tdiss,H 2 °f S as a t the edge of this region. The main 
sink has a LW luminosity of 10 48 s _1 , while the average H2 
column density in the neutral region is ~ 3 x 10 19 cm" 2 . We 
then g et an expression similar to equation 9 of lJohnson et al.l 
(2007), but differing by a factor of ten because we consider 
a less luminous star: 

From this we find that gas at a distance of 0.6 pc will in- 
deed have tdiss,H 2 — 5000 yr, the time of our simulation. 
Disk shielding leads to a morphology of the neutral region 
quite similar to that of the ionized region, as can be seen in 
Figure [5] which shows the expansion of the neutral region as 
molecular inflow becomes more confined to the disk plane. 



3.5 Protostellar Mass Growth 

The sink growth in both cases is very similar for the first 
~ 200 yr, up to when the main sinks reach 12 Mq. This is 
simila r to some of the simulations presented in lSmith et al.l 
l|201ll ) in which feedback did not cause a significant decline 
in sink growth until it grew to > 10 Mq. Afterwards, we 
find the deviation in sink accretion history to be significant, 
leading to a final mass of 27 Mq in the 'no-feedback' case 
and 20 Mq in the 'with-feedback' case. 

In our 'no-feedback' case, the average accretion rate is 
5.4 x 10" 3 M yr" 1 (Fig. [TO - For the first 600 yr, M» evolves 
with time as i aC c 4 . After the ejection of the secondary sink, 
the growth rate declines significantly to M* oc f^cc 3 (see Fig- 
ure [TO}. Even after the growth of the disk and envelope is 
halted, tidal torques are able to continuously funnel addi- 
tional mass onto the sink. 

With feedback, the average accretion rate over the en- 
tire simulation drops to 4.2 x 10 _3 Mq yr -1 . The sink grows 
more slowly as M* oc iacc 6 for the first 500 years, and then 
M t oc t ac ° 9 afterwards. To describe the mass growth history 
in a different way, we also find that the average accretion 
rate declines from 3 x 10 -2 Mq yr -1 over the first 500 years 
to 7 x 10 -4 Mq yr -1 over the remainder of the simulation. 
Given this strong protostellar feedback onto the disk, we 
can devise a simple analytical estimate for the final mass 
attainable by Pop III stars as follows: 

M* = M*tteed, (24) 

where, £f cc a is the timescale over which Pop III accretion will 
be essentially shut off. Given that Pop III stars must accrete 
through a disk, and assuming for simplicity the thin-disk 
approximation, 

M* = 3ttEi/, (25) 

where v is the effective viscosity within the disk. We can 
approximate this using 

v = ctcjn, (26) 

wher e a is the disk viscosity parameter (|Shakura fc Sunvaevl 
1 19731) . For gravitational torques within strongly self- 
gravitating disks, a oc (t CO oif2) -1 and typically ranges from 
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Figure 9. Evolution of various H II region properties over time. Solid line is taken from the simulation, while the dashed line is the 
prediction of the self-similar Shu champagne flow solution. Dotted line also shows the I-front evolution found from the test simulation 
in which the ray-tracing was updated at each timestep instead of every five timesteps. (a): Total mass of the H II region, (b): Radial 
extent of the H II region, taken as the average distance between the sink and the ionized particles, (c): Average density of the ionized 
particles, (d): Average temperature of the ionized particles (solid line). Also shown is the average abundance of neutral hydrogen, /hi 
(dash-dotted line). Note how the H II region fluctuates on small timescales while the long-term evolution generally follows the predicted 
analytical solution. 



0.1-1 (e.g. lLodato fc Ricdl200ot ). In cases like the central re- 
gions of our disk closest to the main sink, where t coo \ ~ Q,^ 1 
(see Section 3.2), the value for a should be closer to 1. 
We now have the expression 

M* = 37rEc 2 ^t fccd . (27) 

If we estimate that the feedback begins once the pressure 
wave of warm neutral gas reaches the scale of the disk ra- 
dius, which will occur over the sound-crossing time i SO und, 
we then have t lccd ~ t sound ~ 1000 AU / 7 km s _1 ~ 700 
yr. Using E = 140 g cm" 2 , a = 1, = 3 X 10" 11 s" 1 , 
and Cs = 2 km s _1 , we have a final mass of M, = 20 Mq, 
well approximating the final protostellar mass reached in our 
'with-feedback' case. 

Given the rapid evolution the disk in both cases at later 
times, let us compare our simulations with analytical esti- 
mates of how the sink accretion rate should decline. In the 
'no-feedback' case, at 1000 yr the disk mass has dropped 
to slightly more than half of its value at the time of initial 



sink formation, and the mass continues to decline further 
to one-fourth the intial value. If we also account for an in- 
crease in disk size due to angular momentum conservation, 
then E decreases by nearly an order of magnitude. The disk 
temperature has also dropped from 1000 to 500 K by the 
end of the simulation (see Figure [2} , leading to a decrease 
in Cs. From this we approximate that the sink accretion rate 
should again drop by over an order of magnitude after 1000 
years, which indeed occurs for the 'no-feedback' case. 

In the 'with-feedback' case, by 3000 years of sink accre- 
tion the mass of gas in the disk has declined from 40 Mq 
to ~ 20 Mq, a factor of two. E thus decreases by the same 
amount. The disk temperature also decreases from 1000 to 
600 K. Considering equations (25) and (26) along with the 
further decline of E due to expansion of the disk, this im- 
plies that the accretion rate should decrease by an order of 
magnitude, from 2 x 10 -2 Mq yr -1 over the first 1000 years 
to ~ 2 x 10 -3 Mq yr -1 afterwards. While the main sink ac- 
cretes almost an order of magnitude more slowly than this 
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after 1000 yr (~ 4 x 10~ 4 M Q yr" 1 ), during this same time 
the secondary sink accretes at 10 -3 Mq yr -1 (Fig.JTTJ. The 
total mass of the sinks thus grows at nearly the expected 
rate, but the main sink undergoes 'fragmentation-induced 
starvation' as the secondary sink intercepts the majority of 
the disk inflow, particularly between 1000 and 2000 yr (Pe- 
ters et al. 2010a). After 2000 yr, however, the distance be- 
tween the sinks increases from ~ 200 AU to greater than 
~ 400 AU, and the secondary sink no longer intercepts the 
inflowing mass. In fact, between 2500 and 3200 yr, the op- 
posite happens and mass flow onto the secondary sink is 
temporarily intercepted by the main sink. 

Also after 2000 yr, radiative feedback begins to shut 
off mass flow onto the disk, and mass flow onto both of 
the sinks gradually declines as well. We can extrapolate the 
mass growth of the main sink to later times by using the 
powerlaw fit to the sink's mass after 2000 yr, A/* oc t^;^?. 
KH contraction will likely end between roug hly 10 4 and 10 5 
yr, and at 10 5 yr our fit yields a mass of ~ 30 Mq for the 
largest star. For the secondary sink, a fit of M, oc t®£? to 
the late-time growth also yields a mass of ~ 30 Mq at 10 5 
yr, implying the stellar system may evolve to an equal-mass 
binary. We again point out that we underestimate the Pop 
III mass, however, due to the early use of ZAMS values 
in our protostellar model and incomplete disk shielding on 
small scales. In comparison, the same extrapolation for the 
'no- feedback' case gives an asymptotic mass of ~ 35 Mq, 
indicating that N-body dynamics can also play a role in 
reducing the final Pop III mass. 

We here remind the reader that our sink algorithm al- 
lowed for sinks to merge. The stars represented by the sinks 
may have instead become a tight binary, and such a bi- 
nary may have been disrupted by later close encounters with 
other stars. Without sink merging, the final mass reached 
by the main sink could have differed (see also discussion in 
iGreif et al.ll201ll ). and the total stellar content might not 
have been dominated by a massive binary. However, our 
main sink gained only a small amount of its total mass, 
~ 4 M©, through mergers, and the second-largest sink also 
gained the majority of its mass through gas accretion. In 
the opposite extreme of no sink mergers, a likely outcome 
may have still been a stellar system dominated by a massive 
binary, but also including several much smaller sinks. De- 
termining the actual stellar merger rate of stars that come 
into close vicinity will be left for future work. 



4 SUMMARY AND CONCLUSIONS 

We have performed cosmological simulations of the build-up 
of a Pop III star to determine how LW and ionizing radia- 
tive feedback influences the mass growth of the star and 
the fragmentation of primordial gas. We find that radiative 
feedback will not prevent fragmentation, but will lower the 
final mass attainable by a Pop III star. When accounting 
for feedback, we estimate a Pop III mass of 30 Mq by 10 
yr, once it has reached the ZAMS. This is a lower limit, 
however, due to details in our modeling of the protostel- 
lar evolution and inability to fully account for inner-disk 
shielding. Inclusion of feedback also led to a massive binary 
syste m instead o f a hig her-order multiple like that seen in, 
e.g., IStacv et alj (|201C ) and our 'no- feedback' case, since 
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Figure 10. Effect of radiative feedback on protostellar accretion. 
Black solid line shows mass growth with no radiative feedback, 
while black dotted line shows the 'with-feedback' case. The blue 
solid line is a double powerlaw fit to the sink growth rate for the 
'no- feedback' case, and the red dotted line is a double powerlaw fit 
for 'with-feedback' case. The dashed line shows the growth of the 
second-most massive sink in the 'with-feedback' case. Radiative 
feedback along with fragmentation-induced starvation leads to a 
divergence in the accretion histories in less than 1000 yr, and in 
the 'with-feedback' case the main sink does not grow beyond ~ 
20 Mq in the time of the simulation. 



the feedback quenched disk growth a nd fra gmentation early 
on. In agreement with ISmith et af] (|201ll ). we furthermore 
find that stellar N-body dynamics can also play a significant 
role in the growth of a Pop III star through stellar ejections 
and disk scattering. The final masses reached in our simu- 
lation agree wel l with the recent two-dimensional study by 
iHosokawa et alj |201lf) . which also examined the effects of 
radiative feedback on Pop III mass growth and found that 
accretion was halted at 40 Mq. 

It is interesting to compare our results to that of re- 
cent work by Peters et al. (2010b). They similarly examine 
ionizing and non-ionizing radiative feedback on massive star 
formation, though they study the case of present-day star 
formation, and their initial configuration was different from 
ours in that they began with a 1000 Mq rotating molecular 
cloud core. They find H n regions which fluctuate in size and 
shape as gas flows onto the stars, and that the final mass 
of the largest stars is set by 'fragmentation-induced star- 
vation,' a process in which the smaller stars accrete mass 
flowing through the disk before it is able to reach the most 
massive star. This is in contrast to models in which the final 
stellar mass is set once ionizing rad iation shuts off the disk 
accretion (e.g. lMcKee~T an 2008). In our 'with-feedback' 
case we find that fragmentation-induced starvation occurs 
between ~ 1000 and 2000 yr, while afterwards radiative 
feedback does lead to a further decline in the sink accre- 
tion rate and the disk mass. Soon after the I-front breaks 
out from the sink, the second largest sink accretes much 
more quickly than the first, intercepting a large portion of 
infalling mass that otherwise would be accreted by the main 
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Figure 11. Accretion rate of various sinks in the simulations. 
Red dotted line, shown in both the top and bottom panels, rep- 
resents the main sink for the 'with-feedback' case. Blue solid line 
(top panel) represents the main sink of the 'no-feedback' case. 
Black dashed line (bottom panel) is for the second-largest sink of 
the 'with-feedback' case. Note the generally larger accretion rate 
of the 'no-feedback' case as compared to the accretion rate under 
feedback. After 1000 yr, the accretion rates of the two largest sinks 
in the 'with-feedback' case generally alternate in terms of which 
has the larger amplitude. The growth of the main sink is especially 
low between 1000 and 2000 yr, a result of fragmentation-induced 
starvation. 



sink. Indeed, the final estimated mass of the main sink is 
similar to that found in Peters et al. (2010a), where they 
found the combination of disk fragmentation and radiative 
feedback led to a most-massive star of 25 M©. However, in 
our case radiative feedback grows in importance after 2000 
yr as it eventually slows mass flow onto eit her sink. A sim- 
ilar st udy of current-day star formation by iKrumholz et all 
(2009) found that a prestellar core would similarly collapse 
into a disk that would host a multiple system of massive 
stars, even under the effects of radiation pressure. However, 



they followed smaller average accretion rates over a longer 
period of time (50,000 yr) and found that gravitational and 
Rayleigh- Taylor instabilities would continue to feed mass 
onto the disk and stars. 

We also note that, despite their very similar numer- 
ical setups, the disk evolution and accretion rate of our 
'no- feedback' case differs somewhat from that described in 
IStacv et all l|2010h . Several distinctions between the simu- 
lations, however, explain this. The high-d ensity cooling and 
chemistry is updated from that used in IStacv et al.l |2010| 
(see Section 2.2). We also use an adaptive softening length 
in stead of a single s oftening length for all gas particles as 
in lStacv et al.l (|2010l ). and our criteria for sink accretion are 
slightly more stringent. The main contribution to the dif- 
ference, however, is likely the stochastic natur e of the sink 
partic le dynamics. While no sink was ejected in lStacv et al.l 
(2010), the sink ejection and subsequent rapid velocity of the 
main sink in our 'no-feedback' case altered the disk struc- 
ture, and the main sink would likely have grown to a higher 
mass otherwise. Nevertheless, the final sink masses in both 
simulations were still the same to within a factor of two. 

The radiative feedbac k seen here is much st ronger than 
the analytical prediction of lMcKee fc Tanl (|2008l ). even when 
considering the combined accretion rate of both sinks. They 
found that a Pop III star could grow to over 100 M© through 
disk accretion, as disk shadowing allowed mass to flow onto 
the star even while the polar regions became ionized. Our 
lack of resolution prevents this disk shadowing from being 
fully modeled on sub-sink scales, and the ionizing photon 
emission emanating from the sink edge is likely overesti- 
mated along some lines-of-sight. This especially applies to 
angles just above and below the disk that quickly became 
ionized in our 'with-feedback case', when in reality these 
regions are unlikely to become ionized until sometime later, 
only after the gas tha t is polar the disk becomes ionized first. 
iMcKee fc Tanl (|2008l ) furthermore assumed disk axisymme- 
try, which does not describe the disk in either of our test 
cases. Nevertheless, our shielding prescription does indeed 
keep the disk from becoming ionized in our 'with-feedback' 
case. However, the I-front does not expand in a perfectly uni- 
form fashion along the polar directions, as the disk rotates, 
and the position of the main sink within the disk varies as 
it orbits its companion sink. Thus, different angles will be 
shielded at different times. Once gas along a certain direction 
has been ionized, it may recombine at later times, but LW 
radiation prevents most of this gas from cooling back down 
to below a few thousand Kelvin. LW radiation, combined 
with recombination of ionized gas and gravitational heating 
from the sink, therefore leads to a bubble of hot neutral gas 
that continues to expand in all directions, except within a 
few degrees of the disk plane. Mass flow onto the disk and 
sinks is then greatly reduced. Thus, while our results under- 
estimate the effect of shielding, it still highlights how non- 
axisymmetry will enhance the effects of radiative feedback, 
and the true physical case likely lies somewh ere in between 
our 'w ith-feedback' case and the prediction of lMcKee fc Tanl 
(2008). In a similar vein, non-axisymmetry can also promote 
further disk fragmentation, and this in turn can result in N- 
body dynamics that may provide another means of reducing 
Pop III accretion rates. 

The fragmentation of primordial gas and growth of 
Pop III stars has recently been modeled from cosmological 
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initial conditions with resoluti on reaching nearly protostel- 
lar scales (Clark et al. 2011b. [Greif et alJlioilf ). However, 
though our simulation is less highly resolved, it explores a 
different regime of Pop III growth. The aforementioned stud- 
ies could not follow the mass accretion for more than 100- 
1000 yr, before the protostars had grown beyond 10 Mq , and 
they did not follow the growth of the I- front. Our work thus 
affords a first look, starting from cosmological initial condi- 
tions, at how Pop III growth will continue beyond 1000 yr, 
after the formation of the I- front. 

In both the 'with-feedback' and 'no-feedback' cases, the 
total stellar mass at 1000 yr was similar to the total stellar 
mass found for th e corresponding time in the simulations of 
iGreif et all |201ll ), which had comparatively more refined 
resolution, utilized sinks with smaller sizes of < 1 AU, and 
found fragmentation on scales < 50 AU. Higher resolution 
may reveal fragmentation o n sub-sink s c ales i n our simu- 
lations as well. Nevertheless, I Greif et al.l (|201ll ) found that 
the highest mass of any individual sink at 1000 yr was al- 
ready ~10 Mq, large enough to emit substantial radiation. 
Our maximum sink masses at this time were slightly larger 
but still similar, 17 and 22 Mq. Thus, while our strong feed- 
back yields a lower limit on the total stellar mass of the 
system, the mass of any given sink is an upper limit to 
the mass of an individual protostar within that sink. Radia- 
tive feedback possibly suppresses fragmentation on sub-sink 
scales, as seen within st udies of cold gas hos ting present- 



20071). Primordial 



day star formation (e.g. IKrumholz et al. 
gas is much warmer, however. and lSmith et ail |201ll ) found 
only mild suppression of primordial gas fragmentation under 
non-ionizing feedback from ~ 10 Mq stars. It remains to be 
seen whether primordial gas would undergo significantly less 
sub-sink fragmentation under ionizing feedback from more 
massive stars. 

These remaining issues show that more computation- 
ally expensive and highly resolved simulations will be neces- 
sary in the future to determine the true nature of sub-sink 
fragmentation. Future cosmological simulations will eventu- 
ally bridge this gap by both resolving protostellar scales and 
modeling ionizing radiation, but such a calculation pushes 
current computational limits. 

Even with these caveats, however, our results are largely 
cons istent with observa t ions o f massive stars in the Galaxy 
(e.g. iMason et al.lll998l . 120091 ). showing that about 40% of 
O stars have companions in the visual binary regime, cor- 
responding t o typical distances o f ~ 1000 AU (see also 
discussion in IKrumholz et al] |2009). Comp utational studies 
of pr esent-day massive star-fomation (e.g. IKrumholz et al.l 
120091 ) have also found disk structures dominated by a mas- 
sive binary. Since similar disk accretion processes are at work 
in the Pop III regime, our finding of a stellar system domi- 
nated by a massive binary is unsurprising. 

If radiative feedback were typically able to prevent Pop 
III stars from growing to more than a few tens of solar 
masses, this would have several effects on their observational 
signatures. PISNe would be less frequent, as this requires 
a star to grow to greater than 140 Mq. Old, metal-poor 
stars within the Milky Way halo and nearby dwarf galaxies 
may preserve the nucleosynthetic pattern of the first SNe, so 
this may help to explain the lack of PISNe chemical signa- 



tures found in these nearby stars (e.g-jchristlieb et al|20"02 
iBeers fc Christl ieb 2005; Frc bel et al.ll2005l ; lTumlinsonl200e 



but see lKarlsson et al"1l2008r i. Instead, Pop III stars may end 
their lives through core-collapse SNe or direct collapse to 
BHs. For sufficient stellar rotation rates , the possibility o f 
Pop III collapsar GRBs also remains fe.g. lstacv et~aT]|201ll ). 
The feedback of Pop III stars on their neighboring metal-free 
minihaloes would also be altered, though the details of how 
the mass and formation rate of such 'Pop III. 2' stars would 
be affected remains to be determined by future simulations. 
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